Numerical Solution of the ID Schrodinger Equation: Bloch Wavefunctions 
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In this article we discuss a procedure to solve the one dimensional (ID) Schrodinger Equation for 
a periodic potential, which may be well suited to teach band structure theory. The procedure is 
conceptually very simple, so that it may be used to teach band theory at the undergraduate level; 
at the same time the point of view is practical, so that the students may experiment computing 
band gaps, and other features of band structure. Another advantage of the procedure lies in that it 
does not use specific symmetry properties of the potential, so that the results are generally valid. 
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I. INTRODUCTION 

Most traditional textbooks dealing with Solid State 
Theory at an elementary level discuss the fundamentals 
of band structure calculations, including in most cases 
only a theoretical description of the essential features of 
band structure, usually with a discussion of the Kronnig- 
Penney modeli. Of course, this is of great pedagogi- 
cal value, but presenting only this viewpoint is unsat- 
isfactory, given the great amount of computer resources 
available today, and the relative ease with which one may 
take advantage of the added insight that may be obtained 
from simple numerical calculations. 

Of course, this fact has not been neglected in scientific 
literature; many articles have been devoted to the nu- 
merical solution of the Schrodinger equation in this and 
other journals 2 ^ 4 ! 5 ! 6 ! 7 ! 8 ! 9 ! 10 ! 11 ! 12 . The references cited 
here intend only to be a (biased) sample of the large 
existing literature on the subject of numerically solving 
the Schrodinger equation for one dimensional potentials. 
Most of the articles seem to concentrate on the problem of 
finding bound states^Si^iiiii^, while a relatively minor 
number treat the problem of finding the energy spectra 
for a periodic potentialSiiiSiiil, which is the subject of 
this article. 

The band theory of Bloch electrons in a periodic po- 
tential is basic for an understanding of solid state physics. 
The usual description is given in most texts in solid state 
theory(see, for example Aschcroft and Mermin's classic 
textbooki, which we follow here for the one dimensional 
case. 

Consider an electron (mass m) moving in a periodic 
potential U(x + a) = U(x), in which the potential U(x) 
is written as a superposition of potential barriers v(x), 
centered the points x = ±(n 4- l/2)a: 



U(x)= J2 



v(x ~ na). 



(1) 



The band structure of the one dimensional solid is ex- 
pressed in terms of the properties of an electron scat- 
tering from a single potential barrier v(x). Therefore, 



they (i) write the wavefunction for an electron scatter- 
ing from the left, with energy E = h 2 K 2 /2m, in terms 
of reflection and transmission amplitudes r and t; these 
wavefunction become, for x > a 



e lKx + re~ lKx 
te %Kx , x > a, 



, x < 



(2) 



while, for the electron coming from the right-hand side, 
the wavefunction is 



te 



-iKx 



—iKx 



x <0 

iKx 



■ re 



x > a. 



(3) 



Observe that these results are valid only for a poten- 
tial that is even with respect to x = a/2, otherwise, one 
would have to introduce different coefficients r' and t' in 
the last equation. We assume, of course, that the rest of 
the solutions, i.e., the part corresponding to the region 
< x < a has been found by some procedure. Clearly, 
since VP; and ^ r are two independent solutions of the 
Schrodinger equation corresponding to the same energy, 
the full solution of the Schrodinger equation correspond- 
ing to the periodic solid will be expressed as a the linear 
combination (now for < x < a) 



*(x) = A^i(x) + B^ r (x), 0<x<a. 



(4) 



According to Bloch's Theorem, the wavefunction ^(x) 
satisfies the relation 



^{x + a) = e tka ^{x), 



(5) 



for suitable values of k. The same relations holds for the 
derivative of ^(x), namely 



(6) 



2 



Imposing the conditions above on the wavefunction ^(x), 
gives a relation that may be used to obtain the energy 
vs. wavevector relation E = E(k) = h 2 K 2 /2m, 



cos(ka) 



if - r A 
2t 



jiKa 



l 

2t' 



-iKa 



(7) 



It is then shown that the energy is a periodic function 
of the Bloch wavenumber k, with period G = 2ir/a (the 
reciprocal lattice vector). 



Ek+G — Ek 



(8) 



Let us concentrate on the solution of the equation for 
the Bloch wavefunction u(x). One way to do this may be 
to expand the wavefunction $ in a Fourier series, doing 
the same for the potential, namely 



*(z) = Y. c ^ qx 
<? 

U(x) = U G e lGx 



e~ lGx U(x)dx, 



in which the Uq are Fourier coefficients, and the sum 
over G goes over reciprocal lattice vectors. 

The equation for the coefficients c q is customarily writ- 
ten as follows 

(J^(k - Gf - E^j c k _ G + J2 U G '-Gc k - G > = 0. (9) 

G' 

From a numerical point of view, this equation may be 
slow in converging to an accurate solution, one may need 
to include a large number of terms in the resulting matrix 
equation. However, this method has been successfully 
used to calculate band structures of real solids (3D). 

In the ID case, we have other methods at our disposal, 
which may be better suited to the problem. Unfortu- 
nately, most of the methods available for the ID case 
cannot be taken over to the 3D case, and our proposal is 
not an exception to this (alas!). 



II. NUMERICAL SOLUTION 

Now let us discuss the particular situation that arises 
in ID. A peculiarity of our treatment is that we really 
have no need for the auxiliary function u(x), as it will be 
seen presently. 

Now, to actually find numerical solutions of the 
Schrodinger equation, and the corresponding energies, 
we proceed as follows. For each value of the energy E, 
considered here as a parameter, we find two independent 



numerical solutions, which we call C(x) and S(x). For ex- 
ample, for a fundamental period with endpoints at x = 
and x = a, these functions may be chosen to satisfy 



C(o/2) 
S(a/2) 



1 C'(a/2)=0 
S'(a/2) = 1 



(10) 
(11) 



To numerically obtain the functions C and S we divide 
the integration region in two regions, < x < a/2, and 
a/2 < x < a. We start the integration from the center 
( x — a/2 towards the boundaries of the two regions, 
imposing at the start the continuity of the function ( C 
os S). After this is done, we compute the normalization 
integral, N — L \^(x)\ 2 dx and redefine the function, di- 
viding by \/N . Observe that the wavefunction defined in 
this way, be it C or S, does not satisfy the full bound- 
ary conditions: it is a real function, and it may not even 
be periodic, the boundary conditions are satisfied only 
by the full wavefunction ty(x). For the numerical inte- 
gration, I have used the Numero\*i£ procedure, since it 
is very accurate and simple to use, however, any reliable 
numerical method may be used just as easily, such as 
Runge-Kutta, with or without adaptive step sizeAi 

The functions C(x) and S(x) are actually dependent 
upon the energy E, i.e., C(x) = Ce(x) and S(x) — 
Se(x), but to make the notation less cumbersome, we 
drop the parameter E henceforth, unless it is necessary 
to keep it. The full wavefunction must be expressed as 
a the linear combination of the functions C and S, with 
possibly complex coefficients A and B, the only way to 
satisfy the boundary conditions from our purely real ba- 
sis functions C and S, 



= A ■ C{x) + B ■ S(x). 
The boundary conditions are now 

*(o) = e 4fca *(0) 
vEr'(a) = e 2fea *'(0), 



(12) 



(13) 



which may be written as equations for the coefficients A 
and B, 



A C A + A S B = 
A' C A + A' S B = 0, 



(14) 



in which the coefficients are may be written as (using 
9 = ka/2) 



Ac 
A s 

A' 



e w C(a) - e- t9 C{0) 
e t6 S(a) - e- t6 S(0) 



jo 



S'(a) - 



• e 

e~ 



e C'(0) 
9 S'(0) 



(15) 
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The equation for the energies is then given by the van- 
ishing of the determinant of the coefficients, 

G(E) - A C A' S - A S A' C = (16) 

When, as in our case, the potential is symmetric with re- 
spect to the midpoint of the period, the functions C and 
5 have definite parity about that point, i.e. C(a) = C(0) 
and 5(a) = —5(0). Also, the derivatives have the oppos- 
ing parity, C'(a) = -C"(0) and 5'(a) = 5'(0), therefore, 

A c = 2isin(0)C(O) 

A s = -2cos(6>)5(0) 

A' c = -2cos(fl)C"(0) 

A^ = 2isin(0)S"(O) (17) 

The eigenvalue equation becomes 

G(E) = sin 2 ((9)C(0)5'(0) + cos 2 (6»)5(0)C"(0) = (18) 

These equations are convenient for numerical work, 
since the functions C and 5 are easily computed numeri- 
cally, and the left hand side of the equation may be plot- 
ted, as a function of E, the zeros being the eigenstates 
that we are looking for. 

In the general (non symmetric) case, we can write 

A c = [C(a)-C(0)}cos(8)+i[C(a) + C(0)}sm(6) 

A s = [5(a) - 5(0)] cos(0) + i[S(a) + 5(0)] sin(0) 

A' c = [C'{a)-C'(0)]cos(6)+i[C'{a) + C'(0)]sm(6) 

A' s = [S'(a) - 5'(0)] cos(0) + i[S'(a) + 5'(0)] cos($9) 

Now, let 





= C(a) 


-C(0) 




= C(a) 


+ C(0) 


5 ( _) 


= 5(a) 


-5(0) 


%) 


= 5(a) 


+ 5(0) 




= C'(a) 


-C'(0) 




= C\a) 


+ C"(0) 




= S'(a) 


- 5'(0) 


s {+) 


= S'(a) 


+ 5'(0). 



We obtain, for the real and imaginary parts of the 
eigenvalue equation: 

X[G{E)] = (c ( _ ) 5 ( '_ ) -5 ( _ ) C ( '_ ) )cos 2 (0) 

+ (5 (+) C ( '_ ) -C (+) 5 ( '_ ) )sin 2 (0) 

%[G(E)] = (c { _)S[ +) +C (+) 5(_ ) - 5 ( _)C( +) - 5 (+) C ( '_ ) 
x sin(0) cos(0) 



First, concentrate upon computing the terms of the 
3(G), they become 

C ( _)5| +) + C (+) 5 ( '_ ) = 2[C(a)5'(a) - C(0)S'(0)] 
5 ( _)C( +) + 5 (+) C ( '_ ) - 2[5(a)C'(a) - 5(0)C'(0)] 

Putting all this together, we find 

%(G(E)) = 2sin(0)cos(6») (W[C,S](a) - W[<7, 5](0)) , 

(20) 

where W[C, S\(<f> = C(<f>)S'((j>) - C"(<?!>)5(0) is the Wron- 
skian of the solutions C and 5; moreover, a well known 
theorem tells us that the Wronskian is a constant, as it 
may easily shown, hence, we have shown that the imagi- 
nary part of G identically vanishes 

$S[G(E)] = 0. (21) 
The eigenvalue equation is therefore, 

G(E) = (c ( _ ) 5 ( '_ ) -5 ( _ ) C ( '_ ) )cos 2 (0) 
+ (5 (+) C ( '_ ) -C (+) 5 ( '_ ) )sin 2 (0) 

III. SIMPLE EXAMPLES 

As indicated previously, a numerical calculation such 
as the one proposed here is able to give us all the energy 
bands that we care to compute (with due consideration of 
the limitations of each integration method). Since com- 
puting wavefunctions is relatively 'cheap', we have devel- 
oped a computer program that does the following: 

1. Divides the first Brillouin zone in N intervals of 
size Afc = 2n/Na. 

2. Prompts the user to indicate the energy range, from 
E min = min(V(x), < x < a) to E max . The 
energy interval is divided into Ne subintervals, of 
size AE = (E max - E min )/N E . 

3. For each fixed value of the wavenumber k — fc , 
it computes the normalized basis functions C and 
5. From these, it computes the function G(Ei), at 
each of the Ne points of the energy interval. These 
values are sent to an array (and stored in a file). 

4. The program seeks the roots of the equation 
G(E) = 0, using the values obtained previously as 
a starting point; the energies are then refined by 
a combination of Newton's method (actually, the 
Secant method) and bisection. 
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5. The resulting roots are stored in a file. Next, 
move over to the next value of k in the /c-grid 
(k = k + Ak), and repeat the previous step. Since 
this is done for each value of k = kj = ko + jAk, 
one obtains the band structure directly from this 
calculation. 

As numerical examples, we show the results obtained 
for three simple potentials, with period 2tt: 

• The Kronig-Penney model (with Vo = 1, width a = 
0.5 centered at x — n) 

V(x) — 0, for \x — 7r| > a 

= Vo for \x-tt\< a, (22) 

• the sinusoidal potential (with Vo = 1 also), 

V{x) = Vq{1 - cos(x))/2, (23) 

• and a triangular potential, 

V(x) — Vq(tt — x)[it, for x < 7r 

= V {x-tt)/tt for x>n. (24) 

These potentials are plotted together in Figure Q we 
have chosen them to have the same maxima, minima, 
and fundamental period < x < 2n. These poten- 
tials have all been thoroughly studied in the literature: 
the Kronig-Penney model has solutions that may be ex- 
pressed in terms of sinusoidal and exponential functions. 
The energy vs. wavenumber relation, E(k), however, is a 
trascendental equation; its solution must be obtained nu- 
merically, which somewhat offsets the comparative ease 
with which one obtains the wavefunctions. As described 
in Stiddard's textbook^ ( in his notation: a square bar- 
rier, of width a, period a + b), the eigenvalue equation 
is 



V(x) Simple test potentials 



- 








V / 


- / 
y 








A 
\ 



0.00 1.80 3J9 5.39 7.18 8.98 10.77 12J7 



FIG. 1: The test potentials used in this work: sinusoidal po- 
tential (1 — cos(x))/2, Kronig-Penney potential centered at 
x = 7r, and a triangular potential. 




FIG. 2: First four energy bands for the three model po- 
tentials. The Kronig-Penney potential is shown in the black 
curve, the sinusoidal potential in the blue curve, and the tri- 
angular potential on the green curve. 



IV. FINAL COMMENTS 



7 2 -/3 2 

sinh(^b) sm((3a)+cosh('-fb) cos(/3a) = cos(k(a+b)). 

2/3 7 

(25) 

The triangular potential may also be solved analiti- 
cally, in terms of Airy functions, but the eigenvalue equa- 
tion becomes very cumbersome to evaluate. A similar 
situation occurs for the sinusoidal potential, which solu- 
tions may be expressed in terms of Mathieu functions. 
In all cases, the analytical calculations needed to obtain 
the eigenvalue equation from the boundary conditions 
become very involved; much more involved that the nu- 
merical solutions described here, and the 'icing on the 
cake', so to speak, is the fact that in the end, all has to 
be evaluated numerically. So, why not do it numerically 
from begining to end?. The results of our calculations 
are shown in Figure [3 



First, we have established a simple numerical proce- 
dure that enables us to compute numerically the wave- 
functions and energy vs. wavevector curve, E(k) for any 
periodic potential. The procedure is very simple, yet it 
appears that it has not been discussed previously, en- 
abling us to study the effect that a parameter change on 
the potential has on the energies E(k), band gaps and 
densities of states g(E). The procedure may be imple- 
mented quite simple using Matlab or Scilab. 
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